#include <iostream>
#include <fstream>
#include "ProbDensStat.hpp"
#include "Dynamics.hpp"
#include "Integration.hpp"
#include "FileValues.hpp"

//using namespace std;
std::string data_dir = "../data/";

int main() {
	Integration integration;

	std::ifstream inFile((data_dir+"physage.in").c_str());
	if (!inFile.good()) {
		std::cerr << "Error when opening input file: " << data_dir << "physage.in";
		exit(1);
	}
	integration.readParameters(inFile);
	inFile.close();

	std::ofstream densFile((data_dir+"dist.dat").c_str());
	integration.setDensFile(densFile);

	std::ofstream statFile((data_dir+"stat.dat").c_str());
	integration.setStatFile(statFile);

	integration.evolution();

	//Store the distribution of the completion time
	/*std::ofstream complTimeDistFile((data_dir+"complTimeDist.dat").c_str());
    ProbDensStat<1> complTime(1);
	integration.completionTime(complTime);
	integration.writeParameters(complTimeDistFile);
	complTime.calcStats();
	complTimeDistFile << "#prob: " << complTime.prob() << "\n";
	complTimeDistFile << "#meanT: " << complTime.meanX() << "\n";
	complTimeDistFile << "#varT: " << complTime.varX() << "\n";
	complTime.saveProbDens(complTimeDistFile,0);*/

	//Scans v and D, saving the mean and the variance of the completion time
	/*std::ofstream compTimeFile((data_dir+"compTime.dat").c_str());
	integration.writeParameters(compTimeFile);
	compTimeFile << "v\tD\tmeanX\tvarX\n";
	ProbDensStat<1> complTime(1);
	for (double v=1; v<=2; v+=1) {
		for (double D=0.01*v; D<=1*v; D*=1.1) {
			std::vector<double> vv(integration.getNStages(),v),
					vD(integration.getNStages(),D);
			integration.setV(vv.begin(),vv.end());
			integration.setD(vD.begin(),vD.end());
			integration.completionTime(complTime);
			complTime.calcStats();
			compTimeFile << v << "\t" << D << "\t" << complTime.meanX() << "\t" << complTime.varX() << std::endl;
		}
		compTimeFile << std::endl;
	}*/

	densFile.close();
	statFile.close();
}
